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ABSTRACT 

The  equatorial  ionosphere  contains  imbedded  bubbles  that  rise  though  a  horizontally  stratified  plasma.  The 
motion  of  the  bubbles  are  affect  by  gravity ,  neutral  winds  or  external  electric  fields  which  produce  electric 
fields  in  the  F-Region  density  perturbations  of  the  bubbles.  Exact  solutions  for  the  electric  potentials  are 
derived  assuming  linear  or  circular  symmetry  to  the  density  structures  imbedded  in  the  background  plasma. 
A  wide  variety  of  analytic  solutions  for  electric  potentials  are  found  for  both  density  cavities  and  density 
enhancements.  An  analytic  description  of  a  rising  bubble  can  be  constructed  by  attaching  a  tail  to  the  top  half 
of  a  circular  hole  to  from  the  electron  density  solution.  The  potential  for  this  plume  structure  is  a  weighted 
sum  of  the  analytic  solutions  for  each  separate  piece.  Using  this  electric  potential,  quasi-analytic  solutions 
for  the  transport  of  the  bubbles  are  derived  using  the  continuity  equation  for  the  plasma  with  production  and 
loss  terms  neglected.  The  analytic  models  of  the  electric  fields  provide  incompressible  motion  that  transports 
the  locations  of  ((plasma  cells  ”  but  does  not  change  the  density  of  the  plasma  in  each  cell.  This  Lagrangean 
approach  employs  a  time  dependent  coordinate  mapping  of  the  undisturbed  layer  grid.  Using  internal 
electric  potentials  of  the  bubbles  and  external  polarizations  of  the  F-layer  as  a  whole,  a  transport  model 
yields  tilted  plasma  plumes  that  move  through  the  F-Region.  This  time-dependent  computer  model  provides 
useful  plasma  densities  in  a  fraction  of  the  time  for  fully  numerical  simulations. 

1.0  INTRODUCTION 

The  F-Region  ionosphere  can  become  unstable  if  a  density  perturbation  becomes  electrically  polarized  by 
external  forces  from  electric  fields,  neutral  winds,  and  gravitational  acceleration.  Near  the  geomagnetic 
equator,  gravity  can  act  on  the  plasma  attached  to  the  nearly  horizontal  magnetic  field  lines  to  produce 
unstable  conditions.  After  sunset  when  the  layer  is  lifted  by  ambient  electric  fields,  the  bottom-side  steepens 
and  plasma  bubbles  are  formed.  These  bubbles  rise  through  the  layer  in  response  to  a  Rayleigh-Taylor  type 
instability.  Also,  winds  or  electric  fields  induce  electric  fields  in  both  density  cavities  and  enhancements  that 
cause  distortions  in  the  density  stmctures.  These  distorted  plasma  stmctures  are  responsible  for  degradation  of 
radio  propagation  that  lead  to  navigation  errors  and  outages,  communications  systems  failures  and  radar 
clutter.  The  modeling  of  ionospheric  bubbles  or  density  enhancements  uses  computer  simulations  the 
calculate  the  effects  of  self  generated  electric  fields  (E)  that  are  driven  by  gravity,  neutral  winds  and  external 
electric  fields.  Numerical  computation  of  the  electric  potentials  requires  the  most  time  and  effort  in  the 
bubble  modeling  process.  The  electric  fields  for  these  simulations  can  be  found  numerically  using  direct  or 
iterative  solvers  of  the  non-separable  potential  equations  that  describe  the  self-generated  electric  fields  (e.g. 
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[1]).  The  computational  time  for  solving  for  the  disturbed  ionosphere  is  often  prohibitive  so  analytic  solutions 
to  both  the  transport  and  potential  equations  are  useful.  Exact  analytic  solutions  for  the  electric  potential  can 
be  used  to  test  the  numerical  algorithms  and  to  determine  errors  produced  by  boundary  conditions  and 
numerical  round-off.  The  analytic  solutions  also  yield  insight  into  the  conditions  for  production  of 
ionospheric  bubbles. 
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Figure  1:  Block  diagrams  of  (a)  numerical  and  (b)  quasi-analytic  algorithms  for  ionospheric  bubble 
modeling.  Both  approaches  use  an  equivalent  set  of  equations  but  apply  different  solution 
techniques  and  different  frames  of  reference.  This  paper  focuses  on  the  solution 
of  the  elliptic  potential  shown  in  red  for  the  (b)  quasi-analytic  approach. 


The  computational  and  analytical  techniques  for  simulations  of  equatorial  bubbles  are  compared  in  Figure  1. 
Typically,  numerical  models  of  equatorial  bubbles  follow  the  procedure  illustrated  by  the  block  diagram  given 
in  Figure  la.  A  stratified  model  of  the  F-layer  is  perturbed  by  a  small  density  disturbance.  Gravity  is  allowed 
to  setup  an  electric  potential  in  this  plasma.  The  electric  potential  is  obtained  with  a  numerical  solution  of  a 
non-separable  elliptic  equation  using  a  direct  solver  such  as  described  in  Appendix  A.  Once  the  electric 
potential  is  obtained,  the  plasma  transported  in  response  to  the  electric  fields  for  a  small  time  step.  A  non- 
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dissipative  flux  corrected  transport  algorithm  [9]  is  then  used  for  incrementally  move  the  plasma  disturbance. 
The  process  is  repeated  with  the  generation  of  a  revised  electric  potential  followed  by  more  incremental 
plasma  transport.  All  of  these  processes  are  numerically  intensive  and  can  require  several  hours  of 
computation. 

A  new  approach  for  the  quasi-analytic  model  of  the  equatorial  bubble  is  proposed  using  the  three  steps  in 
Figure  lb.  First,  the  electric  potential  is  defined  by  an  analytic  function  that  gives  self-consistent  expressions 
for  electron  density  (or  Pedersen  conductivity)  structures  that  are  obtained  in  the  presence  of  background 
electric  fields,  neutral  winds  and  gravity.  This  procedure  is  described  in  the  next  section  of  this  paper. 
Second,  the  plasma  transport  is  determined  using  incompressible  motion  from  the  induced  electric  fields.  The 
plasma  transport  is  derived  with  the  analytic  electric  potential  distorting  the  coordinates  with  out  changing  the 
density  in  each  coordinate  cell.  The  third  step  is  to  adjust  the  parameters  in  the  analytic  models  to  that  the 
analytic  solution  given  in  step  1  matches  the  quasi-analytic  results  from  step  2.  The  application  of  the 
transport  and  normalization  processes  will  be  described  in  the  later  sections  of  this  paper. 


2.0  ANALYTIC  MODELS  FOR  THE  ELECTRIC  POTENTIAL  IN  A  DISTURBED 
IONOSPHERE 

The  equatorial  ionosphere  is  commonly  thought  of  as  a  uniform  layer  with  the  occasional  imbedded  structure 
or  bubble.  The  modeling  of  ionospheric  bubbles  uses  computer  simulations  that  calculate  the  effects  of  self 
generated  electric  fields  (E)  that  are  driven  by  gravity,  neutral  winds  and  external  electric  fields.  The 
equations  for  these  simulations  can  be  found  in  a  number  of  papers  including  [2].  For  the  analytic  solutions 
considered  here,  the  background  ionosphere  will  be  uniform  in  the  horizontal,  x-  and  z-directions.  The 
ambient  magnetic  field,  B,  is  aligned  with  the  z-axis.  The  altitude  variations  of  the  undisturbed  ionosphere 
will  be  represented  by  the  function  ne0(y)  where  y  is  the  vertical  coordinate. 

The  layer  becomes  distorted  when  a  small  perturbation  grows  as  electric  fields  provide  incompressible 
perpendicular  motion  at  F-region  altitudes.  These  internal  electric  fields  move  plasma  across  magnetic  field 
lines  with  the  velocity 

ExB  VOxB 


where  v  is  the  velocity  perpendicular  to  the  ambient  magnetic  field  B  and  the  electric  field  E  =  -VO  can  be 
represented  as  the  gradient  of  a  scalar  electrostatic  potential  O.  Other  perpendicular  components  of  velocity 
driven  by  pressure  gradients,  neutral  winds,  and  gravity  can  be  neglected  because  the  plasma  in  the  F-region 
ionosphere  is  magnetized.  This  means  that  the  electron  and  ion  gyro  frequencies  are  much  larger  than  the 
corresponding  collision  frequencies  with  the  background  neutral  gas. 

An  analytic  approach  is  derived  to  solve  for  the  electric  potential  for  localized  electron  density  perturbations 
driven  by  external  vector  fields  Kelley  [6]  gives  the  electric  current  in  the  F-region  as 


J  =  [E  +  E0  +  (U  +  x  B]  cp 

V;„ 


(2) 


RTO-M  P-IST-056 


18-3 


UNCLASSIFIED/UNLIMITED 


UNCLASSIFIED/UNLIMITED 


Quasi-Analytic  Models  for  Density  Bubbles 
and  Plasma  Clouds  in  the  Equatorial  Ionosphere 


ORGANIZATION 


11  V  ■ 

where  Gp  = - is  the  Pedersen  conductivity,  ne(x,  y,  z)  is  the  electron  density,  Vin  is  the  ion  neutral 

b0  n, 

collisions  frequency,  E0  is  the  external  electric  field  perpendicular  to  B,  U  =  Uxx  +  Uyy  gives  the  normal 

components  of  the  neutral  wind  vector,  g  =  -g0y  is  the  gravitational  acceleration,  B  =  B0z  is  the  magnetic 

field  vector  and  the  inertial  terms  in  the  momentum  equation  have  been  neglected.  Assume  that  the  normal 
electric  fields  are  constant  along  the  magnetic  field  in  the  z-direction  and  that  there  is  no  z-component  of 
current  J.  For  specific  gravitational  accelerations,  neutral  winds  and  external  electric  fields,  the  potential 
equation  is  found  from  (2)  using 

dJ 

-^  +  V±.J  =  0  (3) 

8z 

Substitution  of  (2)  into  (3)  and  integration  along  the  z-direction  leaves  an  equation  for  the  potential  in  the 
perpendicular  (_L)  x-  and  y-  directions. 

V,  •  (Sp  V±0>)  =  IpV±2(D  +  V±<J>  •  V±Zp  =V±-J[E0+(U  +  -^-)xB]  cpdz  (4) 

^in 

where  Zp  =  J  Gpdzis  the  field-line-integrated  Pedersen  conductivity  and  E  =  — V±®  is  the  induced  electric 

g 

field.  For  simplicity,  the  driving  fields  E0,U,  and —  are  assumed  constant  in  space  and  time,  then  the 

vin 

potential  equation  simplifies  to 


V>  =  {[E0  +  (U  +  — )xB]-  V±0}  •  V±  ln(£  )  =  {Ex  -  V.O}  •  V±  ln(Sp)  (5) 

vin 


g 

where  the  equivalent  electric  field  vector  defined  by  ET  =  E0  +  (U  H - )  x  B  .  Given  a  spatial  distribution 

Vin 

for  the  Pedersen  conductivity  (or  electron  density),  the  potential  is  usually  obtained  numerically  from  the  non- 
separable  elliptic  equation  (5).  Often  iterative  solvers  requiring  relatively  long  solution  times  or  direct  solvers 
requiring  large  memories  are  required  to  compute  this  solution. 

A  computational  alternate  approach  assumes  that  the  potential  is  given  and  (5)  is  used  to  find  the  associated 
electron  density.  For  this  solution,  only  Pedersen  currents  in  the  horizontal,  x-direction  will  be  considered  so 

[E0x  +  (U  -  — )B0]  x  =  ETxx  (6) 

V 

in 

where  ETx  represents  the  equivalent  driving  fields  from  a  static  electric  field  in  the  positive  x-direction,  a 
neutral  wind  in  the  positive  y-direction  and  gravitational  acceleration  along  the  positive  y-axis.  The  sign  of 
ETx  is  negative  for  the  downward  acceleration  of  gravity.  In  our  notation,  the  growth  rate  for  the  Rayleigh- 
Taylor  instability  is  y  =  (-ETx/B0)/LN  where  LN  is  the  scale  length  of  the  gradient  on  the  bottom  side  of  the 
ionosphere  [8]. 
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In  normalized  Cartesian  coordinates,  (4)  becomes 


32®  <3Log(2  )  ao  a20  3Log(S  )  ao  3Log(2  ) 


ax2 


-+ 


-+ 


ax  ax  ay 


+  - 


ay  ay 


8x 


=  0 


where  often  gravitational  forcing  is  the  sole  contribution  to  ETx,  ®  = 


® 


Ejx  r0 


(7) 


is  the  dimensionless,  normalized 


potential,  (x,  y)  = 


f  \ 

A  1 

Vro’roy 


are  the  normalized  coordinates,  and  r0  is  a  constant  scale  factor  for  all  distances. 


Note  that  (7)  has  many  self-similar  elements.  Multiplying  the  x  -  and  y  -  coordinates  as  well  as  ®  by  a 
constant  scale  factor  (i.e.,  r0)  does  not  change  the  equation.  Multiplying  Sp  by  a  constant  C0  also  yields  a 

solution.  Consequently,  if  ®(x,  y)  and  Sp(x,y)  satisfies  (7)  then  so  do  the  pairs  of  functions  r0  ®( — , — ) 


ro  ro 


x  y 


and  C0  Sp  (  — ,  —  ) .  Normalized  coordinates  ( x ,  y )  will  be  used  to  simplify  the  notation  for  the  analytic 
solutions. 


ro  ro 


The  existence  of  analytic  solutions  for  (7)  was  discovered  by  examining  numerical  solutions.  A  numerical 
algorithm  for  non-separable  elliptic  equations  similar  to  (7)  was  written  using  a  block  tri-diagonal  solver  of 
the  algebraic  equation  derived  from  finite  difference  approximations  to  the  partial  derivatives  (Appendix  A). 

When  a  circularly  symmetric  function,  £p(r)  where  ?  =  ^/x2  +  y2  ,  was  used  for  the  Pedersen  conductivity  it 

was  found  that  the  integral  of  the  resulting  electric  potential  along  x  was  also  circularly  symmetric.  In 
mathematical  terms, 


2p  (? )  =>  J  <D(x,  y)  dx  =  F,  (r) 


(8) 


This  immediately  shows  that  the  form  of  the  potential  is  the  x  coordinate  multiplied  by  a  circularly 
symmetric  function  since 


F^r)  F^r)  dr  F^r)  x 

1>(x,y)  =  =  F2(r)  x 

ox  dv  ox  or  r 

Functions  of  single  variables  such  as  r  can  be  solved  analytically. 


(9) 


General  analytic  solutions  to  (7)  can  be  derived  by  making  simplifying  assumptions  about  the  form  of  ®  and 
Fog(Sp).  Numerical  simulations  for  symmetric  perturbations  in  x  with  x-directed  fields  given  by  (6)  yield 
electric  potentials  that  are  odd  functions  of  x  with  the  form 


0(x,  y)  =  (a0  +  aj  x)  G[q(x,  y)] 


(10) 


where  a0  and  ai  are  constants,  q(x, y)  =  ^x2+sy2  is  a  single  variable  representing  an  elliptical  perturbation 
for  the  potential  and  “s”  determines  the  polarization  of  the  coordinate  ellipse. 
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The  Pedersen  conductivity  (or  electron  density)  takes  the  form  of  an  elliptically  shaped  perturbation 
modulating  the  background  conductivity. 

Log[SpS)  (x,y)]  =  Lp[q(x,y),s]  +  log[2p0(y)]  (11) 


where  Lp(q,s)  is  the  natural  log  of  the  conductivity/density  perturbation  and  £p0(y)  = 


f  e  ne0(y,  z)  vin 

J  B  Q, 


r0dz 


is  the  integrated  conductivity  associated  with  the  horizontally  stratified  plasma  layer.  The  last  term  of  the  left 
side  of  (7)  is  x-direction  gradient  of  the  log-Pedersen  conductivity  which  drives  the  solution  for  the  potential. 
Consequently,  the  x  variation  through  the  function  q  is  required  to  obtain  useful  solutions  for  the  potential. 


Substitution  of  (10)  and  (11)  into  (7)  and  solving  for  the  derivative  of  Lp(q)  yields  the  equation 


q’G'(q){-  23|S 


Lp(q,s)  =  — 


a0  +ajX 


-  +  s  +  sy 


51ogPp0(y)]i 


cy 


}  +  q3G"(q)  +  y2  (s  -  l)s[qG"(q)  -  G'(q)] 


xq  ta'G(qj  11 + qJ0’(q) + y2(s - i)sqG'(q)] 

a0 


(12) 


where  the  prime  (')  denotes  the  derivative  with  respect  to  q.  If  the  functions  of  y  vanish,  (12)  may  be 
integrated  directly.  The  y 2  terms  vanish  only  if  s  =  0,  or  1 . 

General  solutions  for  an  extended  vertical  plume  imposed  on  a  horizontally  stratified  ionosphere  are 
considered  first.  For  s  =  0  the  potential  has  no  variations  in  the  altitude  coordinate  (y )  and  the  solutions  for 
(10)  and  (1 1)  are  given  by 


q  =  x 


O(0)  (x)  =  (a0 
4(q,0)  =  - 


+  atx)  G(q) 

2  aj  x  G'(q)  +  q  (a0  +  ajx)  G"(q) 
x  [  a!  G(q)  - 1]  +  q  (a0  +  atx)  G'(q) 


y(0) 


(x,y)  =  C0Sp0(y) 


_ -1 _ 

(a0  +  dix  x)  G'(x)  +3lx  G(x)  -  1 


(13) 


where  C0  is  a  constant  chosen  to  give  the  background  conductivity  as  x  — »  oo .  Exact  solutions  for  pairs  of 
Pedersen  conductivity  and  electric  potentials  from  (13)  are  easily  found.  Table  I  gives  several  examples  these 
pairs. 

Exact  solutions  for  pairs  of  Pedersen  conductivity  and  electric  potentials  from  (13)  are  easily  found.  The  last 
equation  in  (13)  can  be  inverted  to  give 


G(x)  =  - 

(ao+a^) 


/[* 


S<°>(q>-']dq 


(14) 
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where  Ep0)(x)  = 


C0  s„(y) 


is  the  normalized  1 -dimensional  density  cross  section  in  the  horizontal  direction. 


From  the  second  equation  in  (13),  the  electric  potential  is  found  as 


6(0)(x)=  J[l  -  £p0) (q)"1  ] dq  (15) 

0 


Table  I.  Electron  Density  Disturbances  and  Companion  1-D  Potential  Functions 


Pedersen  Conductivity  Function,  £p0)(x,y) 

Electric  Potential,  0(x,  y) 

£p0(y)exp(bx2) 

(a0+aj  x)exp(-bx2) 

exp(bx2)  +  2b(a0  +aj  x)x-aj 

£p0(y)(i+|x|b)2 

al  x 

(1+  |  x  |b)2  +  aj(b  - 1)  |  x  b  -al 

1+  I  X  |b 

Spo(y)  cosh[b  x] 

(a0  +a!x)sech[b  x] 

cosh[b  x]  -  aj  +  b(a0  +  ajx)  sinh[b  x] 

An  example  of  the  second  function  in  Table  I  is  illustrated  in  Figure  2  with  the  parameters  ai  =  -0.5  and  b  =  3. 
The  conductivity  trough  (Figure  2a)  centered  along  the  y-axis  has  a  ridges  that  increase  in  amplitude  as  the 
parameter  “af’  is  increased.  If  ai  >  0,  the  trough  is  replaced  by  a  Pedersen  conductivity  enhancement  as  the 
sense  of  the  potential  (Figure  2b)  is  reversed.  The  parameter  “b”  simultaneously  controls  the  steepness  of  the 
walls  on  the  conductivity  irregularity  and  the  spatial  decay  of  the  potential  function.  With  separation  of 
variables,  even  more  general  solutions  to  (7)  can  be  found  in  Cartesian  coordinates.  The  1-D  conductivity 
structure  and  electric  potential  represent  the  “tail”  portion  of  an  ionospheric  bubble. 
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Figure  2.  One-Dimensional  Pedersen  conductivity  associated  with  the  1-D  analytic  electric 
potential  using  the  parameters  ai=  -0.5  and  b  =  3.  The  topology  of  the  solution 
remains  unaffected  by  the  choice  of  model  parameters. 


General  solutions  for  circular  holes  in  an  ionosphere  with  simple  vertical  structure  are  considered  next.  For  s 
=1  the  density  disturbance  is  circularly  symmetric  with  radius  r  around  the  (x ,  y)  origin.  The  solutions  from 
(12)  require  a0  =  0  and  ai  =  1  with  the  result 


r  =  A/x2  +y2 
<l>(1)(x,y)  =  x  G[r] 


{3  +  y 


dlog|Xpo(y)] 


Lp(?,l)  =  ■ 


dy 


}G'(f)  +  f  G"(f) 


G(r)  +  f  G'(?)-l 


(16) 


To  eliminate  any  dependence  on  of  y  in  (16),  the  background  Pederson  conductivity  takes  the  functional  form 


Ep0(y)  =  C0ym  where  m  and  C0  are  a  constants.  With  this  substitution,  (16)  becomes 


4(?,i)  =  - 


(3  +  m)G'(r)  +  r  Gff(r) 
G(f)  +  rG'(f)-l 


(17) 


which  is  identical  to  (13)  with  a0  =  0,  ai  =  1  and  m  =  -1.  The  corresponding  formula  for  the  spatial  variation 
of  Pedersen  conductivity  is 


2(1) 


(£,y)  = 


c0  yn 


G(r)  +  r  G'(r)  - 1 


exp 


-(l  +  m)J 


G'(f) 


G(r)  +  r  G'(r)  - 1 


-dr 


(18) 
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Using  separation  of  variables,  a  wider  range  of  solutions  for  (7)  in  cylindrical  geometry  can  be  derived  but 
( 1 8)  embodies  the  useful  solutions  for  equatorial  bubble  modeling. 

The  one-dimensional  (s=0)  and  two-dimensional,  circularly-symmetric  (s=l)  expressions  are  similar.  With  a0 

cl 

=  0  and  ai  =  1,  the  rational  polynomial  function  of  the  form  G(q)  = - ^  used  in  (10),  and  the 


i+q 


corresponding  potentials  from  (13)  and  (14)  are 


d>l0)(x) 


x  a 

iW 


and  0(1)(x,y)  = 


x  a 

l  +  (x2  +y2)b/2 


(19a) 

(19b) 


where  a  and  b  are  constants.  Analytic  solutions  can  be  obtained  from  (13)  through  (16). 
Derivatives  of  the  density  perturbations  become 


L'(q,0)  = 

L'(r,l)  = 


a  b  qb_1[l  +  b  +  (1  -  b)qb] 

(1+  qb)[a-a(b-l)qb  -  (1  +  qb)2] 
abrbl[2  +  b  +  m  +  (2-b  +  m)rb] 
(1+  r b )  [a-a(b- 1  )r b  -  ( 1  +  r b  )2  ] 


where  q  =  |x| 

=4 


where  r 


x2  +y2 


(20a) 

(20b) 


The  corresponding  Pederson  conductivity  expressions  are 

(l+|x|b)2 


2^(x,y)  =  ^(y) 


1  -  a  -  Bj  |x|b+|x| 


|2b 


2(p1)(x,y)=  exp 


I  a(l  +  m) 


2tan" 


B,  -  2rfc 


•  sign(a)  n 


Li  J 


C0  ym  (l  +  ?b)2 
1  -  a  -  Bjrb  +  r2b 


(21a) 

(21b) 


where  A,  =  yj—  a[4b  +  a(b  - 1  )2 ]  and  B,  =  -a(b  - 1)  - 2  and  r  =  yjx2  +  y2  .  The  constants  of  integration 

are  chosen  to  yield  the  background  density  at  large  distances  where  x  — »  oo .  The  physically  acceptable 
solutions  have  b  >  0. 


The  analytic  solutions  give  insight  into  the  relationships  between  localized  density  perturbations  and  the 
associated  electric  potentials.  Two  types  of  conductivity  structures,  cavities  and  enhancements,  are  described 


by  (21a  and  21b).  In  the  parameter  range 


-4b 


^Min 


>  a  >  0  ,  the  plasma  structure  is  a  cavity  centered  at 


(b-iy 

x  =  0.  These  limits  are  found  by  solving  for  Ai  =  0.  As  parameter  “a”  approaches  the  value  of  “amin”  the  sides 
of  the  density  cavity  becomes  steeper.  With  a  =  amin,  the  wall  of  the  cavity  is  located  at  radius  is 


r  = 


b+i 

b-l 


.  For  b  >  m  +  2,  the  cavity  has  a  ridge  located  at  r  = 


b  +  m  +  2 


-m 
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With  0  >  a  >  1,  a  conductivity  enhancement  is  found  at  the  origin.  This  enhancement  can  represent  the 
increased  Pedersen  conductivity  produced  by  an  artificial  ion  cloud  from  Barium  or  similar  material  released 
in  the  sunlit  ionosphere.  As  a  — »  1 ,  the  sides  of  the  density  enhancement  become  steeper.  Solutions  exist  for 
all  values  of  b  >  0  but  if  b  >  1,  the  potential  vanishes  at  large  distances.  If  a  >  1,  then  the  solutions  in  (21a  and 
21b)  become  complex  and  are  not  physically  possible.  The  maximum  upward  velocity  for  the  potential 


<J>(1)(x,y)  in  (17b)  is  V  0  = 


1  d®(x,y,t)  aE 


Tx 


yo 


B0  dx  B0 


If  a  >  1,  then  the  conductivity  would  move  with  a 


velocity  larger  than  the  ETx/B0  velocity  of  the  driving  force,  which  is  not  possible.  For  instance  if  only  a 
vertical  wind  Uy  is  considered  in  the  equation  (6)  for  ETx,  then  the  upward  velocity  is  VyO  =  a  Uy.  An 
unreasonable  value  of  a  >  1  would  permit  the  conductivity  enhancement  to  rise  faster  than  the  neutral  wind 
driver. 


The  restrictions  on  the  ranges  for  the  potential  amplitude  “a”  indicate  that  not  all  electric  potentials  correspond 
to  a  physical  density  or  Pedersen  conductivity  structure.  For  a  given  force  on  the  plasma  from  external 
electric  fields,  neutral  winds  or  gravity,  the  induced  potential  is  determined  by  the  gradients  on  the  wall  of 
density  cavity  or  enhancement.  These  gradients  are  physically  limited  by  infinite  steepness  and  the  amplitude 
of  the  potential  is  a  maximum  at  this  limit.  Thus,  for  the  solution  to  (4),  a  given  physical  density  structure  will 
always  correspond  to  a  potential  function.  The  magnitude  of  a  potential  function  can  be  increased  to  the  point 
that  there  is  no  corresponding  plasma  density  function. 

The  one  dimensional  expression  (21a)  can  represent  a  horizontal  density  modulation  that  uniformly  changes 
the  background  density  of  a  stratified  ionosphere.  These  may  be  produced  by  horizontally  traveling  acoustic- 
gravity  waves  can  act  as  seeds  for  equatorial  bubbles.  The  elongated  shapes  of  these  modulations  are 
illustrated  by  the  example  in  Figure  2.  The  elongations  can  be  found  in  nature  as  the  extensions  of  an 
ionospheric  plume  below  its  top.  The  horizontal  electric  field  vectors  calculated  as  gradients  of  the  potential 
yield  vertical  plasma  transport.  This  transport  is  normal  to  the  density  gradients  and,  consequently,  no  net 
change  in  the  densities  is  produced.  The  electric  fields  near  the  top  of  the  bubble  are  the  primary  drivers  for 
plasma  transport.  The  cylindrical  solutions  describe  the  head  portion  of  the  bubble.  Combining  the  head  and 
tail  portions  of  the  analytic  bubbles  is  discussed  in  Section  3. 
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a  =  0.5 


Figure  3.  Analytic  results  for  density  cavities  and  enhancements  in  a  uniform  background 
(i.e.,  m=0).  The  densities  and  potentials  are  computed  using  (a)  b  =  2  ,  a  =  0.99  aMin  =  -5.92,  and 
(b)  b  =89,  a  =  0.5  aMin  =  -0.48,  (c)  b  =  4,  a  =  0.5.  The  changes  in  the  parameters  yield  either 
(a)  a  cavity  with  steep  sides,  (b)  a  ridge  around  the  cavity  or  (c)  a  peaked  enhancement. 
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The  two-dimensional  solutions  to  the  potential  equation  are  more  useful  than  the  one-dimensional  solutions. 
The  expression  for  y)  in  (21b)  describes  a  plasma  disturbance  with  two-dimensional  structure.  The 

conductivity  (or  electron  density)  vanishes  in  Sp1}(x,y)at  y=0  unless  m=0.  Figure  3  illustrates  three 

examples  of  the  analytic  density  cavities  and  the  associated  electric  potential  for  a  uniform  background  using 
m  =  0.  By  changing  the  parameters  in  the  analytic  model,  a  wide  variety  of  density  structures  is  obtained. 

Pedersen  Conductivity  Normalized  Potential 


a  =  -1.78 

h  -  A 


a  =  -1.78 
b  =  4 


Figure  4.  Density  cavities  imbedded  in  a  linear  and  inverse-quartic  variation  for  the  vertical 
profile  of  the  plasma  conductivity.  Both  structures  yield  the  same  electric  potential. 

The  conductivity  and  the  y-directed  electric  field  go  to  zero  at  the  lower  boundary. 
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The  background  plasma  variation  can  be  approximated  using  nonzero  values  of  m.  The  Pedersen  conductivity 
from  (21b)  vanishes  at  the  y  =  0  boundary  if  m  >  0.  With  SP  =  0  at  the  lower  y=0  boundary,  the  potential 

ad 


equation  (7)  reduces  to 


dy  dy 


=  0  .  To  satisfy  this  condition,  5®/<9y  =  0  because  a  vertical  gradient 


a£p/<9yis  nonzero-positive  at  the  lower  boundary.  This  condition  is  automatically  built  into  the  analytic 
expression  for  the  electric  potential  because  of  the  y 2  symmetry  of  ®(x,  y)  =  x  G(^/x2  +  y2)  . 


Two  quantitatively  different  ionospheres  can  become  polarized  with  the  same  potential  variation.  Figure  4 
illustrates  two  solutions  of  (21b)  with  identical  parameters  for  a  and  b  but  with  (a)  m  =  1  for  a  linear 
background  profile,  and  (b)  m  =  0.25  for  a  forth-root  of  y  profile.  As  seen  by  the  solutions  in  (21a  and  21b),  a 
family  of  Pedersen  conductivity  structures  can  be  associated  with  a  single  electric  potential  (or  field) 
distribution.  This  non-uniqueness  property  can  be  exploited  for  modeling  the  evolution  of  the  density 
structures 


3.0  FORMATION  OF  BUBBLE  STRUCTURES  FROM  ANALYTIC  SOLUTIONS 
FOR  TROUGHS  AND  HOLES 

The  analytic  solutions  derived  in  Section  2.0  are  the  building  blocks  for  analytic  descriptions  of  equatorial 
bubble  structures.  A  single  plasma  bubble  or  the  individual  finger  of  plasma  bifurcation  can  be  decomposed 
into  holes  and  troughs.  By  adjusting  the  parameters  of  the  analytic  descriptions  for  the  holes  and  troughs,  any 
finger  of  a  plasma  irregularity  can  be  approximated.  The  procedure  for  bubble  finger  formation  is  based  on 

noting  that  an  analytic  potential  in  (14)  with  the  form  0(1)(x,y)  =  x  Gx[^jx2  +  y2  ]  is  shown  to  have  a 
density  structure  of  the  form 


Cn 


V  (x,y)  = - -  exp 

G^rJ  +  fGj  (r)-l 

where  r  =  ^x2  +  y2  and  m  =  0. 


G/(?) 


G^f)  +  f  Gj  (?)  - 1 


-dr 


(22) 


Analytic  bubbles  or  fingers  are  formed  by  using  (22)  for  the  top  portion  or  “tip”  of  the  finger  above  some 
altitude  y0  and  by  extending  the  horizontal,  x-axis  cross-section  downward  for  the  plume  portion  of  the  bubble 
for  y  less  than  y0.  The  bubble  is  then  fully  described  by 


SBubble(Xy) 


|  2j)(x,y)fory  >y0 
l2?)(x,y0)fory  >y0 


(23) 


Away  from  the  transition  altitude  y0,  the  corresponding  electric  potential  is  given  by 


^  Bubble  (X.  ?)  = 


0>(1)(x,  y)  =  xGJ^x2  +  y2  ] 
O(0)(x)  =  j[l  -  S(P1)(q,y0)-1]dq 


for  y  »  y0 


for  y  «  y0 


(24) 
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These  equations  satisfy  the  normalized  potential  equation  (7)  at  all  locations  near  y  =  y0  where  the  electric 
potential  from  one  type  solution  penetrates  into  the  region  of  the  other  type  solution. 

The  bubble  potential  can  be  thought  of  as  a  weighted  sum  of  the  analytic  solutions  for  the  potentials  of  the 
trough  and  hole.  If  the  end  of  the  trough  joins  the  center  of  the  hole  (i.e.,  y0  =  0)  then  the  potential  at  that 

SLogP^x,  0)] 


point  is  the  average  of  the  two  potentials.  This  is  proven  by  noting  at  y  =  0  that  - — 

Qy 

therefore - ^Bubble  =  0 .  Without  the  y-derivative,  the  potential  equation  (7)  becomes 

dy 

^  ^Bubble  _|_  dLog(£Bubble )  dOBubble  |  d  O 

Bubble  _  dLog(£Bubble)  _  q 

dx2  dx  dx  dy2  dx 

By  assuming  that  y0  =  0,  the  potential  equations  for  the  hole  and  the  trough  take  similar  forms 

’dO(1)(x,  0) 


=  0  and 


0"<Du,(x,  0)  a2OUJ(x,  0)  5Log[E^(x,  0)] 


dxz 


dyz 


dx 


dx 


-1 


=  0 


(26) 


and 


a2O(0)(x)  02O,O|(x)  5Log[E^(x,  0)] 


ax2 


dy2 


■  +  ■ 


dx 


aow(x) 

dx 


-1 


=  0 


(27) 


Adding  (26)  and  (26)  and  comparing  with  (25)  immediately  yields  that  the  bubble  potential  as  the  average  of 
the  trough  and  the  hole  potentials 


_  <D(1)(x,0)  +  <Dw(x)  _ 


^BubbleO^0)  = 


(0)/ 


xG1[|x|]  + j[l  -  41)(q,y0)-1]dq 


2  2 
In  general,  the  bubble  potential  from  a  finger  structure  (23)  takes  the  form 

^Bubble (x»  y)  =  (I)"'(x,  y)  g(x,  y)  +  d)(0|(x)  [1  -  g(x,  y)] 
where  the  transition  function  g(x,  y)  is  bounded  by  0  and  1  with  the  limits 
f  0  for  y  «  y0 


(28) 


g(x,  y)  -> 


1  for  y  »  y0 


(29) 


and  for  yo  =  0,  g(x,0)  =  14. 
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Figure  5.  Electric  potentials  for  a  plasma  bubble  (a)  formed  by  attached  a  plume  to  a  circular  density  cavity. 
The  numerical  solution  of  the  potential  equation  (b)  is  approximated  by  the  analytic  solutions  of  the  potentials 
for  the  plume  attached  to  the  potential  for  the  cavity  (c).  The  potential  fluctuations  at  the  tip  of  the 
bubble  closely  match  the  analytic  dipole  potential  (d)  associated  with  the  plasma  hole. 


To  illustrate  the  analytic  finger  formation,  rational  polynomial  function 

G,(q)  =  7^-iT  (30) 

1  +  q 

is  used  to  generate  density  and  potential  structures.  The  bubble  will  be  formed  in  a  uniform  background 
plasma  by  letting  m  =  0.  The  transition  from  plume  to  tip  will  occur  at  y  =  y0  using  (21b)  to  give  the  analytic 
result 

SBubble(X>y)=  eXP 

where  q  =  ^/x2  +  y2  for  y  >  y0  and  q  =  ^x2  +  yf2  for  y  <  y0 


2tan_1 


-2q 

A, 


bA 


sign(a)  7i 


b  \2 


(l  +  qB) 


1-a-Bjq  +q 


2b 


(31) 


RTO-M  P-IST-056 


18-15 


UNCLASSIFIED/UNLIMITED 


UNCLASSIFIED/UNLIMITED 


Quasi-Analytic  Models  for  Density  Bubbles 
and  Plasma  Clouds  in  the  Equatorial  Ionosphere 


ORGANIZATION 


The  electric  potential  is  found  from  (24)  with  Ob  bbl  (x,y)  =  -  for  y  >  y0  and  the  numerical 

l  +  (x2+y2)b/2 

X 

integral  of  Jf1  -  2Bubbie(^yo)_1]dq  for  y  <  y0 .  Figure  5  illustrates  the  bubble  structures  for  the 

o 

parameters  a  =  -0.9,  b  =  3,  and  the  normalized  y0  =  -1.  The  hybrid  structure  of  a  trough  attached  to  a  hole 
yields  a  realistic  description  for  a  plasma  bubble  (Figure  5a).  Using  a  numerical  solver  for  the  potential 
equation  (see  Appendix  A),  the  electric  potential  for  the  analytic  bubble  structure  is  computed  and  is  shown  in 
Figure  5b.  The  hybrid  potential  (Figure  5c)  agrees  very  well  with  the  computation  in  all  regions  except  at  the 
interface  y  =  y0  in  (23)  and  (24)  where  there  is  a  discontinuity  in  the  derivative  of  the  model  bubble.  For 
reference,  the  analytic  potential  associated  with  just  the  tip  of  the  bubble  is  given  in  Figure  5d.  It  turns  out 
that  this  bubble  tip  potential  is  all  that  is  required  for  computation  by  the  transport  equations  in  the  plasma. 

The  bubble  potential  drives  the  plasma  transport  to  change  the  density  of  the  plasma.  The  incompressible 
continuity  equation  is  given  by 


+  v- V 


A 


=  o 


(32) 


where  the  plasma  velocity  is  computed  from  the  electric  potential  using 


V^xB 


(33) 


Normalizations  of  the  potential,  derivatives,  velocities,  and  time  are  given  by 

_  cp  —  v  B  ^  — E 

,  V_l  =  r0V1?  v  =  —  0  ,  and  t  =  —  Tx  t  with  the  negative  sign  for  the  electric  field  used 


0  = 


Ejx  ro 


"'Tx 


B0  ro 


because  with  downward  gravity  ETx  <  0.  With  these  normalizations,  the  continuity  and  plasma  drift  velocity 
combined  to  give  the  rate  of  change  for  the  integrated  Pederson  conductivity 


es  V , O x B  -  /-  ~  ~  \  B 


dt 


B 


o 


(34) 


o 


The  gradient  of  the  potential  (i.e.  electric  field)  must  be  orthogonal  to  the  gradient  in  the  Pederson 
conductivity  for  the  densities  to  change.  In  the  trough  portion  of  the  plasma  bubble,  the  density  gradients  are 
aligned  with  the  gradients  in  the  potential  (Figure  2)  and  there  will  be  no  change  in  the  plasma  density  in  this 
region  where  y  <  y0.  Consequently,  the  primary  driver  of  evolution  for  the  plasma  bubbles  is  the  electric 
potential  near  the  tip  of  the  bubble  structures. 
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Figure  6.  Temporal  changes  in  Pedersen  conductivity  by  incompressible  transport  from  the  induced 
electric  fields.  The  (a)  numerical,  (b)  hybrid,  and  (c)  analytical  potentials  from  Figure  5b,  c, 
and  d,  respectively,  are  used  to  compute  the  rate  of  change  for  the  hybrid  bubble 
model  shown  in  Figure  5a.  The  analytic  potential  for  a  plasma  hole  fit  at 
the  tip  of  the  bubble  yields  a  good  match  to  the  full  numerical  solution. 

Figure  6  provides  a  comparison  for  the  conductivity  temporal  changes  of  (34)  for  the  model  parameters  used 
in  Figure  5  for  a  numerical  potential  (Appendix  A),  hybrid  analytic  potential  model  (24),  and  potential  for  the 
circular  portion  of  the  bubble  (19b).  All  three  descriptions  of  the  potential  yield  a  density  reduction  at  the  top 
of  the  bubble  because  as  the  bubble  rises,  the  low  density  portion  of  the  bubble  is  transported  upward.  The 
density  increase  at  the  bottom  of  the  density  hole  also  represents  a  vertical  rise.  The  transport  of  the  tail 
portion  of  the  bubble  does  not  provide  a  change  in  the  density  or  conductivity.  In  the  case  of  Figure  6,  the 
electric  potential  from  the  circular  part  of  the  bubble  should  be  scaled  up  by  1.2  to  provide  the  correct  plasma 
convection. 

From  this  study,  it  is  concluded  that  it  is  sufficient  to  only  consider  the  electric  potential  associated  with  the 
top  portion  of  the  plasma  bubble  to  compute  its  temporal  evolution.  The  potential  of  the  circular  portion  of 
the  bubble  is  most  easily  found  by  fitting  a  plasma  density  hole  such  as  given  by  (21b)  and  then  using  the  fit 
parameters  to  describe  the  potential  with  the  analytic  representation  of  (19b).  This  procedure  provides  the 
basis  for  the  plasma  transport  algorithm  developed  in  the  sections  4  and  5.  This  discussion  justifies  the  quasi- 
analytic  approach  illustrated  in  Figure  1. 

4.0  QUASI-ANALYTIC  SOLUTIONS  FOR  PLASMA  TRANSPORT  IN  A 
DISTURBED  IONOSPHERE 

Plasma  densities  evolve  by  transport,  compression,  production  and  loss.  These  processes  are  contained  in  the 
continuity  equation 
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— -  +  V  •  (nv)  =  P  -  L  (35) 

ot 

where  n  is  density,  P  is  production  and  L  is  loss.  Numerical  or  analytic  solutions  to  (30)  can  be  based  on 
either  the  Eulerian  or  the  Lagrangean  form  of  the  equations  [7].  These  forms  are  equivalent  except  the 
Lagrangean  form  describes  where  each  bit  of  fluid  cam  from  originally.  Equation  (30)  is  the  Eulerian  from  of 
the  equation  of  continuity.  The  Lagrangean  from  is  derived  for  the  special  case  of  incompressible  flow. 

The  compressibility  of  a  plasma  is  given  by  the  term  V  •  v  .  When  this  term  is  zero,  the  plasma  is  termed 
incompressible.  Using  (1)  the  compressibility  of  the  F-region  plasma  is  found  to  vanish  because 

V  •  v  =  -V  •  (VcD  x  B)/B02  =  [-B  •  (V  x  VO)  +  V(D  •  (V  x  B)]/Bq  =  0  (36) 


where  a  constant  B  is  assumed. 


Expanding  (30)  with  (31)  yields  the  compressionless  form  of  the  continuity  equation 


—  +  v- Vn 

8t 


(  $ 

—  +  v 
V<9t 


\ 


J 


Dn  _ 

n  = =  P-L 

Dt 


(37) 


where  Dn/Dt  is  called  the  total  derivative.  The  total  derivative  of  the  density  moves  with  a  small  volume 
element  (Ax,  Ay,  A z)  in  the  velocity  field.  During  this  process  a  plasma  element  at  (x,  y,  z)  is  mapped  to 
another  location  (x5,  y’,  z’)  in  an  increment  of  time  At .  The  incremental  mapping  equation  is  given  by 

x'  =  x  +  v(x)A^  (38) 


where  the  vectors  have  components  x  =  (x,  y,  z)  and  v  =  (vx,  vy,  vz).  The  derivative  form  of  (33)  simplifies 
(32)  so  that 


dn(x,  t) 
dt 

dx(t)  _ 
dt 


=  P(x,t)-L(x,t) 
v(x,  t) 


(29) 


If  production  and  loss  are  neglected,  then  the  electric  fields  simply  move  volume  elements  of  plasma  in  space 
but  the  density  in  each  element  remains  unchanged.  Equations  (24)  are  the  Lagrangean  form  of  the  continuity 
equation. 

An  analytic  formulation  is  developed  to  describe  the  equatorial  bubbles  in  terms  of  a  mapping  function  that 
distorts  the  ionospheric  layer  according  to  the  second  equation  in  (24).  The  mapping  function  is  usually 
determined  with  a  numerical  simulation  that  calculates  the  electrostatic  E  fields  as  a  function  of  time  and 
space  as  an  ionospheric  bubble  or  irregularity  is  formed  in  the  F-layer.  Substitution  of  these  fields  into  (1) 
yields  the  transport  velocities  and  the  second  equation  in  (24)  can  be  solved  to  provide  the  motion  of  the 
density  coordinates.  The  transformation  of  coordinates  by  this  process  is  given  by 

x(t)  =  M[x(t0),t-t0]  and  x(t0)  =  M_1[x(t),t  - 10]  (30) 
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and,  neglecting  production  and  loss,  the  electron  densities  are  given  by 

n(x,t)  =  n(M_1[x,t-t0],t0)  (31) 


where  the  map  M  1  transforms  the  distorted  coordinates  back  to  the  initial  coordinate  locations. 

Taking  the  magnetic  field  to  be  aligned  with  the  z-direction,  the  electron  density  fluctuations  are  assumed  to 
vary  in  the  2-dimensional  coordinate  system  (x,  y).  The  coordinate  transform  map  is  given  as 


'x(t)- 

Mx(x0,yo,t-to) 

_y(t)_ 

1 

S 

o 

o 

<— K 

1 

t— K 
O 

(32) 


where  y  is  altitude,  x  is  zonal  distance  in  a  flat  earth  system  and  (x0,  yo)  are  the  initial  values  for  these 
coordinates  at  time  t0. 

The  mapping  function  M(x,  y,  t)  must  be  one-to-one  and  invertible  and  be  the  identity  map  where  the  induced 
electric  potential  is  zero  and  the  plasma  densities  are  unchanged.  Substitution  of  (1)  into  (29)  yields 


dx(t) 

dt 


VO  xB 

Bn 


(33) 


Assuming  uniformity  along  B  in  the  z-direction,  the  differential  equations  governing  the  coordinated 
transformations  are 


dx(Xo,y0,t)  _ _ 1  dQ(x, y,t)  ^  dy(x0,y0,t)  _  1  dQ(x,y,t) 


0t 


Bn 


dy 


dt 


Bn 


dx 


(34) 


Differentiating  (34)  by  xO  and  yO  and  using  the  Poisson’s  equation 

a2o  a2o 

V-E  =  V.(-V4.)  =  -(-t+-t) 

5x0  Sy, 

yields  the  equations  for  the  velocities 
d(dx/dt)  _  d(dy/dt)  ,  d(dx/dt)  d(dy/dt)  _ 


dx, 


o  dy0  dy0  — 0 


dxn 


(35) 


(36) 


The  areas  between  curves  of  constant  xO  and  yO  are  preserved  in  the  transformation  and  the  trajectories  of  the 
coordinate  transformation  follow  contours  of  constant  0(x,y).  The  electric  potential  is  setup  by  gravity, 
neutral  winds  and  external  electric  fields. 

Consider  the  coordinate  transformation  provided  by  the  analytic  potential  from  (19b)  in  the  previous  section 

(X/ro)  ax  r0  ETx 


0(x,y)  = 


l  +  [(xAo)  +(y/r0)  ] 


2  -,b/2 


x  ax  ro  ETx 

1  +  ?b 


(37) 
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where  coordinates  x  and  y  are  the  Cartesian  coordinates  normalized  by  a  constant  scale  factor  r0  and 
constants  ax  and  b  define  the  shape  of  the  potential. 

The  plasma  transport  velocities  are 
8x  bxy  r b~2 

W  ~  (i+fb)2 

(38) 

8y  _  1-bx2  r b"2  +  r b 

at-  (1  +  rb)2 

where  t  =  —  ^Tx  t  and  r  =  Jx2  +  y2  . 

r0B0 


The  maximum  upward  velocity  at  the  ( x  =0,  y  =  0)  origin  is  given  by  Vy0  = 


a  E 


Tx 


Bn 


independent  of  the 


potential  shape  parameter  “b”.  Note  that  in  our  example  of  the  downward  gravity  vector  driving  the  transport, 
the  parameter  ETx  from  (6)  is  less  than  zero  and  a  value  of  ax  <  0  is  required  to  yield  an  upward  velocity.  The 


parameter  ax<0  denotes  a  density  cavity  and  consequently  the  center  of  the  cavity  is  expected  to  rise  against 
gravity.  If  the  parameter  “ax”  were  greater  than  zero,  the  center  of  the  density  enhancement  would  fall  as 
expected  under  the  influence  of  gravity.  For  the  rest  of  the  discussion,  only  density  cavities  will  be 
considered. 


The  sides  of  a  cavity  fall  under  the  influence  of  gravity.  The  minimum  downward  velocity  of 

/i  /i  .  i\l/b 


a  ET 

v  =_  x  Tx 

v  wi 


(b-i)2 


yi 


=  -v. 


(b-1)2  •  ,  ,  .  -  J b  +  T 

-  is  found  at  x  =  ±  - 

lb-1 


and  y  =  0.  The  Cartesian  coordinate 

B  4b  yu  4b  U>-i; 

system  is  distorted  by  the  potential  inside  a  conductivity  cavity  so  that  the  center  cells  move  upward  and  the 
side  cells  move  downward. 


The  analytic  model  potential  with  b  =  4  was  previously  illustrated  in  Figure  3b.  The  corresponding  vector 
field  for  the  plasma  velocities  (Figure  6a)  shows  the  central  uplift  of  the  plasma.  As  a  result  of  this  flow,  after 

normalized  time  t  =  1  the  initially  square  cells  become  mapped  according  to  the  results  shown  in  Figure  7b. 
The  horizontal  (red)  and  vertical  (blue)  grid  lines  become  distorted  by  the  vortex  flow  from  the  potential. 
Note  that  the  area  in  each  plasma  cell  remains  constant  during  this  process. 
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Figure  7.  Computed  (a)  plasma  velocities  and  (b)  mapping  of  Cartesian  coordinates  by  the  analytic 
model  for  the  electric  potential  shown  in  Figure  3b.  The  longest  velocity  vector  has  a  magnitude 

a  E. 

- —  and  the  spatial  coordinates  are  normalized  by  the  scale  length  r0.  The  coordinate  map  is 


illustrated  after  the  flow  velocities  flow  have  operated  on  the  plasma  for  time 


t  = 


VO 


-i 


or  when  normalized  time  t  =  1 . 


The  primary  feature  of  ionospheric  bubbles  is  that  they  rise  leaving  an  elongated  cavity  of  reduced  plasma 
density  often  referred  to  as  a  “plume”.  As  the  plume  rises,  it  carries  the  along  electric  potential.  Using  the 


upward  velocity  for  the  potential  in  (32)  ofVy0  = 


av  E 


Tx 


Bn 


,  the  analytic  model  for  the  potential  becomes 


0(x,y,t)  =  <D(x,y-VyOt,O)  = 


x  ax  ETx 


1  + 


/  \ 


vroy 


+ 


y.V'i 

Vro  ro  j 


-ib/2  * 


(39) 


With  this  potential  formula,  the  plasma  drift  velocities  are  found  from  (34)  to  be 
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dx  _  ~ 

ar=x(y-,) 


b[ 

x2+(y-t)2] 

b-2 

2 

t* 

[x2  +  (y- 1)2 

f) 

dy 

dt 


l  +  [(l~b)  x2  +  (y- t)2][x2  + 


(y-<)21 


b-2 

2 


(40) 


where  the  time  has  been  normalized  by  t  =  t  Vy0/r0  and  [x(0),  y(0)]  =  [x0,y0]  are  the  initial  conditions  at 

time  t  =  t0  =  0 .  This  equation  is  the  base  for  the  quasi-analytic  transport  for  production  of  ionospheric 

bubbles.  By  using  this  equation,  it  is  assumed  that  the  electric  potential  shape  is  totally  specified  with  fixed 
parameters  ETx,  a,  b,  and  r0.  The  only  temporal  variation  of  the  potential  is  that  it  rises  with  a  constant  speed 
given  by 


Vx 


y0 


a 


The  potential  in  the  trailing  portion  of  the  plume  is  neglected  because  any  vertical  flow  below  the  bubble  only 
operates  on  horizontal  gradients  and,  consequently,  there  is  no  net  plasma  transport  within  this  region. 


In  the  reference  frame  of  the  rising  potential,  y  =  y  —  t  and  the  differential  equations  become 


The  coordinate  transformation  if  found  by  integrating  (41)  for  (x,y  )  with  an  initial  value  of  (x0,y0)  and 

then  finding  y  =  yp  -  t  .  Further  simplification  is  obtained  by  transforming  to  spherical  coordinates  where 
x  =  rp  Cos0p  and  yp  =  rp  Sin  0P  .  With  this  substitution,  (41)  becomes 


5rp  _  rp  sin  0P 
dt  1  +  ?pb 

30p  _  rpb  l  cos0p(l  +  b  + rp) 

~8f=  + 

where  the  initial  conditions  at  t  =  0  are  rP(0)  =  and  0P(O)  =  tan_1(y0/x0)  . 


(42) 
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The  coordinate  transform  starts  with  an  electric  potential  at  altitude  y  =  t0  =  0  and  lets  this  potential  distort 

the  coordinates  until  time  t,  when  the  potential  is  at  altitude  y  =  \ .  The  starting  and  stopping  times  and 
altitudes  are  critical  in  defining  the  coordinate  distortions.  For  describing  equatorial  bubbles,  these  starting 
altitude  must  be  transformed  to  the  location  where  the  bubbles  starts  to  form  on  the  bottom  side  of  the 

ionosphere  and  the  stopping  altitude  yields  with  the  location  of  the  bubble  at  time  ix  =  t,  r0/Vy0 .  This 
renormalization  is  described  later. 

The  properties  of  this  coordinate  transform  can  be  examined  at  the  center  where  x  =  0  and  (40)  simplifies  to 


—  =  0  with  x(0)  =  x(t)  =  0 
dt 


8y 

dt 


— - - -v  with  y(0)  =  y0 

i+[(y-t)!]= 


(43) 


The  solutions  of  (43)  are  found  from  the  nonlinear  equation 

ly-tr  I  y  |1_b 

y  Sign(y0)  +  =1  y0 1  + '  ,  (44) 

1-b  1-b 


The  limiting  forms  for  the  asymptotic  solutions  to  (44)  are 

y  =  y0for  y0<-l 

y  =  y0  + 1  if  t  <  tc 


y  =  y0  +  tc  if  t  >  t, 

y  =  y0  for  0  <  y0  and  t  <  y0 
y  =  t  for  y0  =  0 
1 


for  - 1  <  y0  <  0  and 


t  -  (~yp) 
c  b-l 


1-b 


y  =  t +■ 


(b-i) 

y  =  t  +  [(b-l)(t-y0)  +  y0 


for  1  <  y0  and  t  =  y0 


,  j_ 

1-b 


for  0  <  y0  and  t  <  y0 


(45) 
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Figure  8.  Mapping  of  the  (a)  y-axis  and  (b)  a  square  coordinate  grid  with  a  rising  potential  that  starts 
at  the  origin.  Using  equation  (34)  with  b  =  4,  the  coordinates  above  the  initial  center  of  the  potential 
are  swept  into  a  narrow  band  that  rises  with  the  upward  bubble  velocity.  The  snapshot  at 

t  =  6  shows  transport  of  the  ionosphere  upward  at  the  center  and  downward  at  the  sides. 


This  maps  of  the  coordinate  distortions  are  illustrated  in  Figure  8  using  the  potential  with  index  b  =  4.  The 
temporal  variations  in  the  normalized  altitude  are  given  in  Figure  8a.  The  spatial  mapping  of  the  grid  at  time 

t  =  6  is  shown  in  Figure  8b.  The  ionosphere  is  nearly  undisturbed  below  the  starting  altitude  ( y  =  0 )  of  the 
rising  electric  potential.  The  ionosphere  is  also  undisturbed  for  times  ( t  <  y0 )  when  the  potential  is  below  the 
undisturbed  ionosphere.  As  the  potential  rises,  it  sweeps  up  all  the  coordinates  and  carries  them  in  a  narrow 
layer  at  just  above  the  center  of  the  potential  function  where  y  =  t  .  This  type  of  coordinate  transformation 
can  be  used  to  simulate  the  rising  equatorial  bubble.  At  time  t,  the  potential  rises  to  be  centered  at  an 
altitude  y1  =  t1. 
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Figure  9.  Coordinate  compression  at  the  top  of  the  plasma  bubble.  For  the  potential  parameters 
b  =  4,  the  normalized  y  -coordinates  are  mapped  to  a  thin  shell  about  0.4  to  0.45 

above  the  center  of  the  rising  potential  function  located  at  y  =  t  . 


The  horizontal  coordinates  near  the  center  of  the  rising  potential  map  to  a  thin  shell  centered  at  the  altitude 
y  =  t .  An  approximation  to  this  map  is  given  by  the  last  equation  in  (45).  The  evolution  of  this  coordinate 
“compression”  is  illustrated  in  Figure  9  for  normalized  times  between  t  =  1  and  20  using  the  parameter  b=4 
for  the  electric  potential.  The  offset  ( y  —  t )  from  the  center  of  the  rising  bubble  increases  monotonically  but 
slowly  as  the  initial  coordinate  covers  a  much  larger  range. 


The  normalized  density  gradient  at  the  edge  of  the  bubble  is 
compression  factor.  Analytically,  this  factor  is  given  by 


aie  _  *9ne  dy0 
dy  dy0  dy 


so 


dy0 

dy 


is  the  coordinate 


dy 


(46) 
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This  factor  is  plotted  in  Figure  10  with  b  =  4  for  a  wide  range  of  times  and  initial  coordinate  altitudes.  The 
compression  factor  increase  with  time  easily  attaining  values  greater  than  10  or  100.  With  this  compression,  a 
bottomside  ionospheric  gradient  becomes  greatly  amplified  as  the  bubble  rises. 


/  ° 

(>^ 


Figure  10.  Density  gradient  enhancement  factor  for  the  coordinate  compression  at  the  top  of  the 
rising  bubble  with  index  b  =  4.  The  compression  becomes  enhanced  with  increases  in  the 

normalized  time  t  .  Outside  the  effect  of  the  potential  this  factor  is  unity.  This  processes 
yields  the  steepened  gradients  at  the  sides  of  the  ionospheric  bubble. 


To  obtain  the  electron  density  or  Pedersen  conductivity  at  time  ( t ),  cell  coordinates  must  be  obtained  for  the 
original  cell  that  was  transported  to  the  current  position.  The  model  ionospheric  bubble  is  therefore  formed  by 
operating  on  a  horizontally  stratified  layer  with  the  inverse  of  the  coordinate  transformation  defined  by  (40). 
This  inverse  may  be  obtained  with  by  (1)  interpolation  of  the  solution  to  (40)  or  (2)  by  reversing  time  for  the 
solution  of  (35).  Figure  8  shows  the  location  of  the  coordinates  at  the  current  time  for  each  normalized 

r  \ 


starting  position  (x0,y0)  = 


2^0.  Yo 

9 


The  interpolation  process  is  hampered  by  the  distortion  of  the 


V  ro  ro) 

coordinate  density  cells.  In  the  mapping  shown  by  Figure  8b,  the  square  cell  with  comers  at  (x0,  y0)  =  (0, 
2),  (0.1,2),  (0.1,  2.1),  and  (0,  2.1)  are  mapped  to  the  elongated  region  with  comers  (0,  6.421),  (0.1075,4.374), 
(0.1112,  4.443),  (0,  6.424)  at  t  =  6.  With  linear  interpolation,  the  point  (x,  y)  =  (0.055,  5.4)  inside  the 
elongated  region  is  determined  to  map  to  the  initial  position  (x0,  y0)  =  {0.0505,  2.0364).  This  is  in  error;  the 
correct  location  for  this  position  is  (x0,  y0)  =  (0.3329,  -0.1784).  Beside  being  prone  to  error,  this  technique 
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requires  interpolation  on  a  non-uniform  mesh  or  numerical  solution  for  the  inverse  of  interpolation  on  an  a 
uniformly  spaced  mesh  of  the  coordinate  locations  at  time  t  =  0 .  For  good  accuracy  and  ease  of 
computation,  interpolation  or  numerical  inverse  solutions  should  be  avoided. 

Since  the  electric  potential  does  not  evolve  with  time,  the  inverse  coordinate  mapping  is  easily  achieved  with 
time  reversal.  Consider  a  cell  with  location  (xi,  yi)  at  time  t  =  ti.  Running  time  backwards  to  time  t  =  0  yields 
its  starting  point  (x0,  yo).  If  the  parameters  vy0  and  b  are  constant,  the  time  reversal  solution  is  most  easily 
obtained  by  replacing  t  with  - 1  in  (35)  through  (37).  This  process  yields  the  map  represented  by  (31). 


The  equations  for  the  inverse  coordinate  map  are 


^=_-o(-o+?)_bk±&±*£j 


b-2 

2 


l  +  [x02 +(y0  + 1)2] 


dt 


1  +  [(1 


b)  Xq  +  (y0  +  t)2  fc  +  (y0  + 1 )2  ]  2 
|l  +  [xg  +  (y0  + 1)2]2  j 


(47) 


The  center  of  the  potential  function  starts  at  altitude  y0  =  t,  and  then  (47)  solved  as  the  center  of  the  potential 
falls  to  an  altitude  y0  =  0 .  For  this  reason,  the  inverse  coordinate  map  equations  are  integrated  from 
t  =  — tj  to  t  =  0 .  The  initial  conditions  for  (47)  are 

x0(-ti)  =  x1  andy0(-t1)  =  y1  (48) 

The  inverse  coordinate  map  is  illustrated  in  Figure  1 1  for  the  initial  grid  and  the  distorted  inverse  grid  at 
several  times.  This  map  is  used  to  determine  the  origin  of  a  coordinate  cell  and  the  initial  electron  density  or 
Pedersen  conductivity  in  that  cell.  As  an  example  of  using  this  inverse  map,  the  mapping  of  specific  point 
(0.055,  5.4)  is  directly  obtained  with  (47)  to  yield  the  correct  initial  location  (0.3329,  -0.1784). 
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Normalized  Distance,  x 

Figure  11.  Inverse  coordinate  map  showing  source  regions  for  density  gradients. 


No  analytic  solution  to  (47)  could  be  found  but  the  numerical  solution  to  the  coupled  ordinary  differential 
equations  is  fast  and  efficient  using  standard  explicit  methods.  The  Lagrangean  inverse-map  is  consequently 
called  quasi-analytic  because  no  finite-difference  approximations  are  applied  to  the  spatial  (x,  y)  coordinates 
but  discrete  steps  in  the  time  dependent  solutions  of  (47)  is  required.  Applying  the  inverse  map  to  any 
analytical  model  of  a  uniform  ionosphere  can  efficiently  provide  the  density  at  any  point  in  space  and  time 
that  can  be  obtained  without  the  use  of  interpolation  or  the  numerical  solution  of  two-dimensional  partial 
differential  equations. 

As  an  example,  the  inverse  coordinate  map  is  applied  to  an  F-layer  described  by  the  following  formula  for  a 
modified  Chapman  layer. 
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Ne(y)  =  Ne0  Exp[l  -  z  -  Exp(z)] 

y-Hp 


Z  = 


Hn 


H0  =  H01  +  [0.5  +  Tan  ( 


H, 


)  /  7T  ]  H. 


02 


(49) 


The  parameters  Ne0,  HP,  H0,  H0i,  H02  and  Hi  control  the  shape  of  the  layer.  The  analytic  simulation  uses  peak 
density  Ne0  =  106  cm'3,  peak  altitude  HP  =  400  km,  bottom-side  scale  height  H0i  =  20  km,  top-side  scale  height 
H02  =  50  km,  and  transition  scale  Hi  =  10  km.  This  simple  layer  model  has  a  steep  bottomside  representative 
of  the  equatorial  ionosphere.  The  conversion  from  normalized  coordinates  is  y  =  y  r0  +  yc0  where  yc0  is  the 
starting  altitude  of  the  electric  potential  at  time  t  =  0. 
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Figure  12.  Simulation  of  a  bubble  formed  in  the  ionosphere  using  a  simple  formulation 
for  the  electric  potential.  The  parameters  for  the  model  are  b=4,  and  ro  =  30  km. 


RTO-MP-IST-056 


18-29 


UNCLASSIFIED/UNLIMITED 


UNCLASSIFIED/UNLIMITED 


Quasi-Analytic  Models  for  Density  Bubbles 
and  Plasma  Clouds  in  the  Equatorial  Ionosphere 


ORGANIZATION 


Applying  the  inverse  coordinate  transformation  (47)  to  the  ionospheric  profile  yields  the  bubble  evolution 
illustrated  in  Figure  12.  The  coordinate  distances  are  determined  using  a  scaling  factor  r0  =  30  km  and  the 
potential  function  rises  through  the  layer  starting  at  yc0  =  370  km  altitude.  The  potential  function  index  is 

arbitrarily  set  to  b  =  4  for  this  example.  The  normalized  time  coordinate  t  =  t  Vy0/r0  =  t  ax  ETx/(B  r0)  is 

used  because  the  parameter  “ax  has  yet  to  be  specified.  The  allowable  values  for  a  were  previously  given  after 

—  4b 

(19)  as  - —  <  a  <  0  .  With  a  larger  value  of  parameter  “ax  ETx”,  the  vertical  velocity  of  the  bubble 

(b  1 ) 

increases  and  the  absolute  time  (t)  in  the  simulation  is  reduced  for  a  fixed  normalized  time.  Figure  12 
illustrates  that  the  analytic  model  using  the  rising  potential  yields  a  quasi-analytic  solution  that  resembles 
numerical  solutions  requiring  much  more  computation  time. 


Zonal  Distance  (km) 


Figure  13.  Effect  of  electric  potential  index  b  on  the  model  ionospheric 
bubble.  The  normalized  time  for  each  solution  is  t  =  3.75  . 


Whereas  the  parameter  “ax”  controls  the  bubble  rise  rate,  the  “b”  determines  the  size  of  the  bubble.  The 
parameter  “b”  may  have  any  value  greater  than  unity.  Larger  values  of  b  yield  electric  potentials  with  larger 
gradients  at  the  edges.  Increasing  b  increases  the  region  affected  by  the  electric  potential.  Figure  13 
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illustrates  the  ionospheric  bubble  for  several  values  of  this  index.  Values  of  the  potential  index  in  the  range  2 
<  b  <  4  seem  to  yield  reasonable  descriptions  of  the  rising  ionospheric  bubble. 


5.0  BUBBLE  TILTS  AND  AMBIENT  SHEAR  FLOW 

Both  observations  and  numerical  computations  [8]  have  demonstrated  that  neutral  winds  can  tilt  the 
ionospheric  plumes  off  vertical.  When  acted  on  by  a  zonal  neutral  wind  Ux,  a  vertical  electric  field  Ey  is 
produce  by  interactions  with  the  background  electron  density  (or  Pedersen  conductivity).  Also,  a  perturbation 
electric  field  ETy  is  produced  by  polarization  of  the  plasma  density  bubble  by  the  neutral  wind.  These  two 
processes  work  simultaneously  to  affect  the  tilt  the  bubbles.  First,  the  background  neutral  wind  induces  large 
scale  plasma  motions  in  the  zonal,  eastward  direction.  The  shear  in  this  large-scale  plasma  drift  will  cause  the 
regions  of  lower  background  Pedersen  conductivity  to  lag  behind  the  regions  of  maximum  Pedersen 
conductivity.  Second,  polarization  of  the  ionospheric  bubbles  by  the  neutral  wind,  which  moves  faster  than 
the  bulk  motion  of  the  plasma,  will  produce  horizontal  motion  that  will  cause  the  density  depletions  to  move 
in  the  opposite  direction  of  the  neutral  wind  relative  to  the  bulk  plasma  drift.  This  is  a  well  known  property  of 
plasma  holes  as  they  respond  to  neutral  winds  [2].  Both  of  these  processes  are  easily  incorporated  in  the 
quasi-analytic  bubble  model. 

Vertical  gradients  in  the  background  density  yield  vertical  gradients  in  the  induced  electric  field  and  vertical 
shears  in  the  horizontal  plasma  drifts  produced  by  these  fields.  This  process  is  captured  in  (2)  assuming  that 
the  vertical  currents  vanish  with  the  result 


Jy  =  [Ey  -  UXB0]  <Tp  =  0  (50) 

The  field  line  is  divided  into  the  F-region  where  the  wind  is  uniform  and  the  E-region  where  the  neutral  wind 
will  assume  to  vanish  [8].  Calling  the  integrated  Pedersen  conductivity  in  these  two  regions 
£f  and  SE  respectively,  the  vertical  electric  field  profile  is  given  by 


Ey(y) 


UxB0SF(y)=  OP,, 

2E  +  2F(y)  dy 


(51) 


where  O0  is  the  polarization  potential  of  the  background  plasma. 
The  resulting  plasma  velocity  is  given  by  (24)  with  the  result 


8x 

~dt 


sF(y)  u  i 

2E+2F(y)  x  f(y/r0)  + 1 


(52) 


where  f(y/r0)  =  f(y)  =  .  (53) 

£F(y/ro) 

This  is  the  horizontal  velocity  of  the  plasma  in  which  the  bubble  is  imbedded.  Usually  the  wind  shear 
variation  is  small  compared  to  the  average  bulk  motion. 
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In  normalized  coordinates,  this  wind  shear  equation  for  the  background  horizontal  motion  becomes 


5x  _  vxs  _  ux  SF(y/r0) 


1 


-VxS(y) 


et  Vy0  yy0SE+IF(y/r0)  1  +  f(y) 

Integration  of  (49)  yields  the  simple  coordinate  map  from  this  large  scale  plasma  motion 


(54) 


x( t )  =  {  VxS  (y)  dt  =  x0  +  t  VxS  (y) 


(55) 


where  x0  is  the  initial  coordinate  at  time  t  =  0  . 


Neutral  Density  (cm  3)  Electron  Density  (106  cm-3) 


Figure  14.  Profiles  from  neutral  atmosphere  and  electron  density  models  used  for 
the  sample  computations  of  electric  polarization  that  tilts  the  equatorial  bubbles. 
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Conductivity  (m2/Ohm)  (m/s) 

Figure  15.  Model  profiles  of  field-line  integrated  Pedersen  conductivity,  uniform  neutral  wind, 
large-scale  horizontal  plasma  drifts,  and  horizontal  wind  relative  to  the  drifting  plasma. 

If  the  E-region  Pedersen  conductivity  is  zero,  then  the  plasma  will  move  horizontally  with  the  neutral  wind 
speed  Ux.  A  finite  coupled  with  vertical  shears  in  the  F-region  Pedersen  conductivity  gives  a  shear 
structure  to  the  horizontal  plasma  motion.  The  tilts  of  the  equatorial  plume  structures  will  be  computed  using 
the  Lagrangean  approach  to  the  plasma  transport  with  both  imbedded  potential  given  by  (39)  and  the  large 
scale  distortions  described  by  (55). 

As  an  example  of  this  wind  induced  shear  in  the  horizontal  plasma  drift,  a  uniform  neutral  with  Ux  =  100  m/s 
was  used  to  polarized  the  plasma  layer  given  by  (49)  with  a  peak  density  of  106  cm"3.  The  model  neutral 
atmosphere  and  electron  density  profiles  illustrated  in  Figure  14  are  used  to  derive  the  equatorial  profile  of  the 
field-line  integrated  Pedersen  conductivity  shown  in  Figure  15.  A  dipole  magnetic  field  model  of  the  form 
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JH^Cose  H^SM 


(56) 


is  used  for  the  magnetic  field  lines  where  H0  =-0.31 1  10'5  Tesla.  Figure  16  illustrates  the  distortion  of  a  square 
mesh  using  the  coordinate  transform  defined  by  (52).  Application  of  this  transform  after  the  transform 
illustrated  in  Figure  8  will  tilt  the  bubble  to  the  left  (west)  side  of  the  simulation.  The  simulation  for  Figure  16 
used  a  fixed  E-region  conductivity  ZE  that  was  one-tenth  the  maximum  value  of  the  F-layer  Pedersen 
conductivity  SF  . 


100  m/s  Neutral  Wind 


Plasma  Drift  (m/s)  Horizontal  Coordinated 

Displacement  (km) 


Figure  16.  Horizontal  coordinate  distortion  by  wind  induced 
plasma  shear  flow  after  500  seconds  of  plasma  motion. 


The  rising  bubble  will  be  caught  in  the  sheared  plasma  flow  to  provide  a  tilt  to  the  bubble.  To  model  this  tilt, 
the  quasi-analytic  transport  model  will  be  modified  assuming  that  the  bubble  follows  a  trajectory  that  is  a 
combination  of  the  vertical  rise  velocity  VyO  and  the  horizontal  shear  velocity  given  by  (52).  The  equations 
for  the  trajectory  of  the  center  (xc,  yc)  of  the  bubble  is 

f  =  V,and^  =  V,s(yc)  (57, 

Assuming  that  the  bubble  starts  to  form  at  the  initial  location  (x,  y)  =  (0,  yc0),  then  the  solution  of  (57) 
becomes 


Yc  =  YcO  +  VyO1  aIld  Xc  =  J  VSx  (YcO  +  VyoO  dt'  (58) 

0 


18-34 


RTO-MP-IST -056 


UNCLASSIFIED/UNLIMITED 


UNCLASSIFIED/UNLIMITED 

Quasi-Analytic  Models  for  Density  Bubbles 
and  Plasma  Clouds  in  the  Equatorial  Ionosphere 


In  normalized  coordinates  this  solution  for  the  center  trajectory  of  the  bubble  becomes 

t 

yc  =  t  and  xc  =  J  VxS(t)  df  (59) 

0 

where  distance  is  normalized  by  r0,  velocities  are  normalized  by  Vy0  and  time  is  normalized  by  r0/Vyo  as 
before. 

The  simple  coordinate  map  (57)  neglects  horizontal  flow  from  internal  polarization  of  the  bubble.  This  flow 
is  the  result  of  vertical  polarization  fields  generated  by  polarization  of  the  bubble  structure  by  the  neutral 
wind.  This  process  has  been  described  in  Section  II  except  that,  along  with  vertical  gravitational  acceleration, 
the  horizontal  neutral  wind  induces  an  electric  potential  that  causes  horizontal  motion  of  the  bubble  center. 
From  (6),  the  vector  electric  field  from  these  forces  is 


E0  +  (U  +  — )  xB  =  -Et 

V;„ 


(60) 


The  horizontal  neutral  wind  in  the  rest  frame  of  the  plasma  bubble  is  Ux  -  Vx0  where  Vx0  is  the  horizontal 
velocity  of  the  center  of  the  bubble.  The  electric  field  associated  with  this  neutral  wind  in  this  frame  is 

Ely  y  =  (Ux  -  Vx0)B0  y  (61) 


by  the  definition  in  (60).  In  the  rest  frame  of  the  plasma,  the  relative  velocity  of  the  electric  potential  at 
altitude  y  =  yc , 


vxR(yc) 


ETy  _  i  ao, 
B0  B0  dy 


(62) 


where  Oi  is  the  potential  from  internal  polarization. 

The  response  of  the  plasma  to  this  drive  field  is  dependent  on  the  physical  structure  of  the  bubble.  For  the 
tilted  bubble  model  with  internal  polarization,  a  single  parameter  scaling  of  the  ambient  drifts  is  used.  The 
horizontal  drift  velocity  of  the  bubble  center  is  assumed  to  have  the  form 

Vx0(y)  =  (1  -  ay)  VxS  =  VxS  +  VxR  with  VxR  =  -  ayVxS  (63) 


where  ay  is  a  constant  analogous  to  ax  used  in  the  previous  section  to  determine  the  bubble  rise  rate  from 
gravity.  The  parameter  ay  determines  the  bubble  velocity  in  a  rest  frame  of  the  background  plasma. 

The  wind  and  plasma  drift  profiles  in  Figure  15  illustrate  this  polarization  effect.  Figure  15b  shows  that,  near 
the  peak,  the  relative  wind  in  the  plasma  rest  frame  is  about  10%  of  the  total  plasma  drift.  The  parameter  ay 
controls  the  relative  velocity  of  the  bubble  in  the  background  plasma.  If  ay  =  0,  the  bubble  drifts  with  the 
background  plasma  as  if  there  were  on  internal  polarization  of  the  bubble.  If  ay  =  1,  the  background  horizontal 
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plasma  drift  cancels  the  wind-induced  relative  motion  produced  by  the  polarization  fields  inside  the  bubble 
giving  Vx0  =  0  and  VxR  =  -VxS.  In  normalized  coordinates  (63)  becomes 

=  tt~  =  0  -  ay)  V,s  =  VxS  +  VxR  .  (64) 

^y0 


Using  the  definition  of 

ay  =  -VxR/VxS  (65) 

where  VxR  is  the  measured  bubble  velocity  relative  to  the  background  plasma  drive  VxS  . 

The  technique  for  bubble  modeling  is  based  on  the  motion  of  the  analytic  electric  potential  along  a  trajectory. 
The  tilted  bubble  rises  along  the  path  defined  by  the  velocity  V0  =  Vx0x  +  Vy0y  .  With  both  background  drift 

and  internal  plasma  motion,  the  dynamics  of  the  center  of  the  plasma  bubble  potential  are  given  by 

-£■  =  V*  ™d  =  VI0  =  (1  -  ar)v,s  (66) 

with  the  initial  conditions  [xc(0),yc(0)]  =  [  0,  yc0  ]  . 

For  the  tilted  bubble  model,  the  internal  horizontal  and  vertical  electric  fields  are  considered  along  with  of  the 
overall  motion  of  the  background  plasma.  With  the  internal  field  assumption,  the  analytical  potential  function 
becomes 


0(x,  y,  t)  =  <D0  (y)  +  ® !  (x  -  xc ,  y  -  yc0  -  Vy0t,  0) 

[x  -  xc  (t)]  Vy0  -  [y  -  yc0  -  Vy0t]  VxR  (yc0  +  Vy0t) 


=  00(y)  +  Bo 


1  + 


c ~  - -V  f 


x-xc(t) 


v  Ao  y 


+ 


y-yco-vy 


b/2 


(67) 


which  is  identical  to  (39)  with  O0,  VxR  and  xc  equal  to  zero.  The  induced  plasma  bubble  velocity  in  the 
horizontal  direction  is 


i  ao 

B0  dy 


V  = _ _ _  =  V  +V 

v  xO  -  vxS  ^  vxR 


(68) 


at  the  position  x  =  xc(t)  and  y  =  yc0  +  Vy0t .  The  added  variables  tilt  the  electric  potential  off  vertical  so  that 
the  head  of  the  bubble  can  flow  against  the  ambient  drift  of  the  background  plasma. 

The  coordinate  shift  xp  =  x  -  xc(t)  and  yp  =  y  -  yc0  -  Vy0t  translates  the  potential  into  the  reference  frame 
of  the  bubble  center  with  the  result 
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®(x,y,t)  =  ®0(y)  +  B0 


XpVyO~yP  VxRCycQ+VyoO 

1  +  (rp/ro)b 


(69) 


where  rp2  =  x2  +  y2  . 

The  potential  function  given  by  (65)  is  substituted  into  the  Lagrangean  transport  equations 


SxO^yoO  _  1  d®(x,y,t)  and  _  1  3®(x,y,t) 


(70) 


dt  B0  dy  dt  B0  dx 

The  resulting  coordinate  mapping  equations  for  the  rising  potential  in  sheared  plasma  flow  is 


dt 


=  vx0(yc)  =  vx0(yc0  +  vy0t) 

b-2  -b 


dx  b  r;-2r-Byp[xpVy0  -  yp  VxR(yc0  +  Vy0t)]  +  [1  +  rpV]  VxR(yc0  +  Vy0t) 


■  +  rpV 


r 


dt  (l- 

5y  =  Vy0[l  +  rpV]-bryVxp[xpYy0  -  yp  VxR(yc0  +  Vy0t)] 

at 


■  +  YxS(y)  (71) 


(1  +  rp  rob) 


incorporating  both  internal  and  external  forces  on  the  bubble. 
In  normalized  coordinates,  the  transport  equations  become 

^  =  vxo(yc)>  yc  =  b 
dt 

dx=  b  (pl>2y[x  -  VsR  (yc)y]  +  [l  +  rpb]  VsR  (yc )  ~ 

[1  +  f]2  xsYj 


at 


(72) 


dy  1  +  r  -  b  x  [x  -  VsR  (yc )  y] 

at 


b  n2 


[l  +  fpD] 


where  once  again  distance  is  normalized  by  ro,  velocity  is  normalized  by  Vyo,  and  time  is  normalized  by  the 
quatntity  r0/Vyo .  The  velocity  functions  are  all  related  to  the  normalized  plasma  shear  by  the  parameter  ay 

withVx0(y)/(l-ay)  =  - VxR(y)/ay  =  VxS(y)  . 

Before  computing  the  coordinate  mapping,  the  trajectory  of  the  center  of  the  potential  function  (xc,  yc)  is 
found  by  solving  the  first  ordinary  differential  equation  in  (71)  or  (72).  For  this  calculation  and  all  the  follow, 
the  shear  velocity  shown  in  Figure  16  was  used  in  (66)  with  a  several  values  of  ay.  Figure  17  illustrates  the 
model  results  for  the  motion  of  bubble  center  as  a  function  of  the  parameter  ay.  With  ay  =  0  the  potential  will 
rise  and  drift  east  reaching  a  maximum  distance  on  the  topside  where  the  wind  induced  drift  vanishes.  As  the 
parameter  ay  is  increased  toward  unity,  the  internal  polarization  of  the  bubble  inhibits  its  horizontal  motion  in 
the  background  plasma  flow. 
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Figure  17.  Curves  of  the  trajectory  for  the  analytic  potential  function  used  to  create  the  equatorial 
bubble.  With  parameter  ay  =  0,  the  internal  polarization  of  the  bubble  is  neglected  and  the  potential 
function  is  transported  by  the  full  action  of  the  plasma  shear.  As  ay  is  increased,  the  zonal  motion 
of  the  center  of  the  bubble  potential  is  reduced.  With  ay  =  1,  the  internal  polarization  at 
the  head  of  the  bubble  completely  cancels  the  plasma  drift  of  the  ambient  layer. 


To  illustrate  the  results  from  these  coordinate-mapping  equations,  the  parameter  ay  is  set  to  one-tenth  so  the 
potential  function  moves  horizontally  with  the  ambient  plasma  flow  at  each  altitude  and  is  slightly  retarded  by 
internal  polarization.  When  ay  =  0.1,  then  VXR(y)  =  VxS/10,  Vx0  =  9  Vxs/10.  The  initial  conditions  for  the 

coordinate  map  are  [xc(0),  x(0),  y(0)]  =  [0,  x0,y0]  at  t  =  0 .  All  of  the  results  are  displayed  in  normalized 
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coordinates.  The  normalization  parameters  r0  =  30  km,  and  vy0  =100  m/s  can  be  used  to  translate  the 
solutions  back  to  physical  space.  The  bubble  moves  with  the  vertical  rise  velocity  from  internal  electric  fields 
set  up  by  gravity.  By  setting  the  internal  polarization  parameter  ay  nearly  to  zero,  the  horizontal  motion  is 
primarily  with  the  background  plasma  but  there  is  a  small  retardation  from  the  internal  fields.  The  shear 
function  was  put  into  normalized  form  from  the  physics  coordinates  using  the  definition  y  =  (y-yc0)/r0 

where  yc0  =  370  km.  With  the  parameter  b  =  4,  the  ordinary  differential  equations  given  in  (67)  are  integrated 
in  time  to  yield  the  quasi-analytic  solutions  for  the  Lagrangean  coordinate  distortions  shown  in  Figure  18. 
The  center  of  the  potential  function  as  derived  from  the  first  equation  in  (72)  is  shown  by  the  green  curve  in 
each  figure.  The  plume  structure  below  the  top  of  the  bubble  becomes  caught  up  in  the  ambient  flow  to  from 
a  backwards  “C”.  The  successive  images  in  Figure  18  are  normalized  shown  for  normalized  times  of  2,  4,  and 
6.  Using  the  normalization  factor  r0/Vyo  =  300  seconds,  the  absolute  times  for  the  tilted  bubble  coordinate 
maps  are  (a)  600,  (b)  1200,  and  (c)  1800  seconds. 
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Figure  18  Coordinate  map  for  distortions  from  a  rising  (b  =  4)  bubble  in  a  sheared  background 
plasma  flow.  The  green  line  shows  the  trajectory  of  the  center  of  the  bubble  with  ay  =  0.1 . 
The  plume  becomes  curved  as  it  is  caught  in  the  ambient  plasma  flow. 


The  coordinate  transform  equations  in  (72)  yield  the  normalized,  destination  coordinates  [x(t),  y(t)]  from  the 
given  initial  coordinates  [x0,y0].  To  determine  a  mapped  density  at  a  given  location,  the  time-inverse 

transformation  to  determine  the  initial  coordinates  for  a  coordinate  cell  that  is  transported  to  a  given  location. 
This  inverse  transform  has  already  been  discussed  in  the  previous  section  for  bubbles  with  out  tilts.  This 
inverse  map  are  obtained  by  replacing  t  with  -t  in  (72).  The  center  of  the  potential  function  starts  at  location 

X0  =  xc(t,)  and  y0  =  \  ■ 
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The  solution  for  the  inverse  coordinate  transform  proceeds  in  two  steps.  First,  the  initial  equation  in  (72)  is 
solved  to  yield  the  function  xc(t)  for  the  horizontal  displacement  of  the  bubble  potential  function.  Next,  the 

system  given  by  (68)  is  integrated  as  the  potential  follows  a  trajectory  [xc(t),  yc(t)]  to  end  at  the  origin  of 
the  normalized  coordinate  system  where  x0  =  0  and  y0  =  0 .  The  curves  in  Figure  17  show  the 
[xc(t),  yc(t)]  trajectories  as  a  function  of  the  internal  polarization  factor  a y.  As  previously  discussed  with 

(47),  the  inverse  coordinate  map  equations  are  integrated  from  t  =  — tx  to  t  =  0 .  The  initial  conditions  for  (72) 
are 


x0(-ti)  =  x1and  y0(-ti)  =  yi 


(73) 
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Figure  19  Computed  examples  of  (a)  inverse  coordinate  map  at  t  =  4  and  (b)  corresponding 
ionospheric  bubble  densities  at  t  =  1200  seconds  for  the  electric  potential  moving  through 
a  sheared  plasma  with  small  internal  polarization  (ay  =  0.1).  The  parameters  for  the 
simulation  are  identical  to  those  used  to  generate  Figure  12  which  can  be  used  for 
comparison  to  illustrate  the  effects  of  the  zonal  wind  on  tilting  the  bubble. 


The  inverse  coordinate  map  for  the  rising  bubble  the  sheared  plasma  flow  is  illustrated  in  Figure  19a  in 
normalized  coordinates  at  the  normalized  time  t  =  4 .  This  coordinate  map  presents  the  origin  of  the  plasma 
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cells  that  have  been  transported  to  the  west  by  the  wind  induced  horizontal  transport  and  the  gravity  induced 
vertical  transport.  The  parameter  ay  is  set  to  zero  so  the  horizontal  transport  by  internal  polarization  of  the 
equatorial  bubble  is  neglected. 

When  the  inverse  coordinate  map  is  applied  to  the  background  plasma  layer,  the  electron  densities  are  to  form 
the  tilted  plume  structure  in  Figure  19b.  The  depletion  at  the  center  of  the  bubble  is  the  result  of 
incompressible  convection  from  the  bottom  to  the  topside  of  the  layer.  In  this  model,  the  reduced  density 
channel  extends  over  100  km  down  though  the  layer  to  the  bottom  of  the  F-region.  The  absolute  spatial 
dimensions  are  derived  using  a  scale  length  r0  =  30  km  and  a  base  height  yc0  =  370  km.  As  with  Figure  17,  the 
time  normalization  factor  is  300  seconds. 


6.0  NORMALIZATION  IN  THE  ANALYTIC  MODELS  OF  THE  DISTURBED 
IONOSPHERE 

In  the  simulations  of  the  previous  three  sections,  the  electric  potential  was  fixed  as  it  rose  through  the  plasma 
layer  forming  bubble  structures.  The  next  step  in  the  quasi-analytic  modeling  is  to  adjust  the  model 
parameters  to  make  the  computed  electron  density  consistent  with  the  electric  potential.  The  appropriate 
values  of  both  parameters  (ax  and  b)  are  obtained  by  comparing  the  analytic  solutions  given  in  (19)  with  the 
quasi-analytic  solutions  obtained  by  compressing  the  F-layer  coordinates  using  (47)  for  an  untilted  bubble  or 
(71)  for  a  plume  with  internal  polarization  and  plasma  shear  flow. 

The  vertical  bubble  motion  comes  from  polarization  of  the  horizontal  density  gradients.  The  horizontal 
Pedersen  conductivity  from  (21b)  though  the  center  of  the  electric  potential  is  given  by 


Vp1)(x,0)=  exp 


a(l  +  m) 

A, 


2tan_1 


^Bj  -  2  |  x 
v  Ai 


ib  X 


J 


-  sign(a)  n 


C0  (l+|x|b)2 
1  -  a  -  Bj  |  x  |b  +x2b 


(74) 


where  Aj  =  y[-  a[4b  +  a(b  - 1  )2  ]  and  B1  =  -a(b  - 1)  -  2  and  |  x  |  =  VT  .  Assume  that  the  electron  density 

at  the  equator  is  directly  proportional  to  the  integrated  Pedersen  conductivity.  In  absolute  coordinates,  the 
analytic  model  for  the  horizontal  electron  density  profile  through  the  center  of  the  bubble  is 


Ne(xp,yc(t))=  exp<! 
(71) 


a(l  +  m) 


A, 


2tan' 


B, 


■2iyR,  n 

A, 


■  sign(a)  7t 


ib  \  2 


C,  (1+  |xp/R,  |  ) 


l-a-B,|xI/R,f  +  |xp/R1 


i2b 


where  xp  =x  — xc(t)  is  the  horizontal  distance  relative  to  the  center  of  the  potential  function  and  {xc(t), 
yc(t)}  is  the  location  of  this  center.  The  corresponding  electric  potential  is 


O(x,y,t)  =  B0 


[x  ~  Xc(t)]  Vy0(t)  +  [y  -  yc(t)]  VxR (yc(t)) 


1  + 


x~xc(t) 

v  R-l(t)  j 


T-Yc(t) 

R,(t) 


2  — |  b(t)/2 


(75) 


where  ax(t),  b(t),  Ri(t),  and  Ci(t)  are  parameters  that  will  be  allowed  to  vary  with  time  as  the  bubble  evolves. 
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The  bubble  rise  velocity  Vy0(t)  = 


ax(Q  Etx 

Bn 


and  the  bubble  retardation  velocity  VxR  =ayVxShave  been 


defined  in  the  previous  sections.  The  electric  potential  follows  a  trajectory  given  by 

rW 

^  =  Vx0(yc)  =  (l-ay)VxS(yc) 

3yc  _v  _  ax(t)  EXx 

at  ”  yoW  ”  Bn 


(76) 


With  this  trajectory,  the  Lagrangean  coordinate  map 


ox(x0,y0,t)  _  1  aO(x,y,t)+v  ^  and  pyi^Yo^)  _  1  3<I>(x,y,t) 


at 


Bn 


dy 


at 


B0  8x 


(77) 


is  again  used  to  determine  the  distortion  of  the  plasma  layer.  The  normalization  procedure  fits  the  function 
given  by  (73)  to  the  densities  determined  using  the  Lagrangean  transport  given  by  (77)  as  applied  to  the 
stratified  plasma  profile. 


Tests  of  the  normalization  procedure  has  demonstrated  that  the  numerical  value  for  the  temporal 
normalization  constant  ax  ETx/(B0  r0)  is  equal  to  approximately  70%  of  the  calculated  growth  rate 


y  =  (-ETx/B0)/Ln  of  the  instability  so  that  the  parameter  ratio  -r0/ax  «LN/0.7.  The  recommended 

procedure  for  providing  reasonable  models  of  equatorial  bubbles  is  to  first  choose  a  scale  length  r0  that 
matches  the  dimension  of  the  “seed”  needed  to  produce  the  bubble.  Second,  select  the  potential  amplitude 


0.7  r0/LN 


using  the  simple  expression  ax 
bottomside  of  the  background  ionosphere. 


where  LN  =  Ne(y)/ 


3Ne(y) 

8y 


is  the  initial  scale  length  of  the 


7.0  CONCLUSIONS 

At  this  point,  all  the  steps  outlined  in  Figure  lb  have  been  completed  and  the  densities  for  the  ionospheric 
bubble  can  be  computed  with  relative  ease.  The  only  numerical  computation  is  a  solution  of  the  ordinary 
differential  equations  (72)  which  are  applied  to  the  model  of  the  unperturbed  ionosphere.  Physics  based 
simulations  of  equatorial  plumes  can  be  computed  with  relatively  high  speed  using  a  simple  form  for  the 
electric  potential  that  has  4  parameters  that  are  adjusted  based  on  numerical  model  fits  to  the  electron  density 
at  each  time  step.  This  fit  procedure  occurs  only  at  the  altitude  of  the  center  of  the  bubble  potential  function. 
The  spatial  and  temporal  scales  for  the  simulations  are  all  normalized  with  a  constant  r0.  The  Lagrangean 
Map  can  be  applied  to  a  wide  variety  of  ionospheres.  This  flexibility  allows  the  generation  of  a  wide  range 
equatorial  bubbles  without  complete  re-computation  of  the  densities.  The  examples  have  illustrated  the 
integrated  Pedersen  conductivities  plotted  at  the  equatorial  plane.  The  densities  along  the  magnetic  field  lines 
may  be  obtained  by  solving  for  one-dimensional  field-aligned  diffusion  of  the  plasma  as  they  are  transported 
by  the  bubble  electric  fields. 

The  formulas  described  here  provide  the  electron  densities  that  can  be  used  for  a  wide  range  of  ionospheric 
applications  including  ray  path  propagation,  diffraction  screen  formation,  radar  and  navigation  error 
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estimation,  and  communications  systems  effects  prediction.  The  major  use  of  this  model  is  to  predict 
observations  from  new  sensors  launched  into  space.  The  results  of  this  model  have  been  predictions  of  the 
signature  of  equatorial  bubbles  on  GPS  occultation  receivers  in  low-earth  orbit  (LEO)  [3]  and  the  phase  and 
amplitude  scintillations  that  would  be  recorded  from  ground-to-space  propagation  from  UHF  beacons  to  LEO 
satellites  [4].  A  secondary  use  of  this  formulation  is  to  provide  incite  into  the  formation  of  equatorial  bubbles. 
The  analytic  formulation  has  shown  that  (1)  a  continuous  range  of  electron  density  structures  will  yield  the 
same  electric  potential  distribution,  and  (2)  the  electric  potential  that  contributes  to  the  evolution  of  the  bubble 
is  concentrated  near  the  head.  Future  research  using  this  technique  will  examine  the  triggering  of  bubbles 
using  initial  distributions  shown  in  Figure  4  and  bifurcations  of  bubbles  where  the  single  potential  function  is 
analytically  split  into  a  pair  of  potentials. 


APPENDIX  A.  NUMERICAL  SOLUTION  OF  THE  POTENTIAL  EQUATION 
USING  A  DIRECT  SOLVER. 


Numerical  solutions  are  required  when  conditions  of  simplified  geometries  or  uniform  flows  yield  a 
complicated,  nonlinear  partial  differential  equation  for  the  electric  potential.  In  Cartesian  coordinates,  the 
non-separable  elliptic  equation  that  describes  the  electric  potential  is  given  by  (7)  and  is  repeated  in  equivalent 
form  here 

*1^05  ^85.®,  (A1) 

p  dx2  dx  dx  p  dy2  dy  dy  dx 


The  conductivity  function  Zp(x,  y)  is  known  and  the  potential  function  0(x,  y)  is  found  as  a  numerical 
solution  to  the  equation  (Al). 

This  equation  is  converted  into  a  set  of  linear  equations  using  the  usual  finite  difference  approximations  to  the 
derivatives  given  by 


a2o  ^-2^+0,^ 

dx2  Ax2 

g2$  ci>,i ,  -2<i>l,+ci>ii-l 

dy2  Ay2 

dx  2Ax 

d®  o,,l+l  -o,,„ 

dy  2Ay 


(A2) 


The  Pedersen  conductivity  is  sampled  in  the  uniform  solution  grid  spaced  by  Ax  and  Ay  to  form  the  array 
variables  2pijwith  (i=l,2,  ...,  M)  and  (j=l,  2,  ...,  N).  To  complete  the  solution,  boundary  conditions  of 

periodic,  fixed/Dirichlet,  derivative/Neumann  or  mixed  form  are  provided. 


RTO-MP-IST-056 


18-43 


UNCLASSIFIED/UNLIMITED 


UNCLASSIFIED/UNLIMITED 


Quasi-Analytic  Models  for  Density  Bubbles 
and  Plasma  Clouds  in  the  Equatorial  Ionosphere 


ORGANIZATION 


The  unknowns  Oi  ■  for  (i=l,2,  . . M)  and  (j=l,  2,  . . N)  are  grouped  into  linear  arrays  given  by 
-  ®mj ]  (A3) 

The  resulting  linear  system  can  be  written  as  an  extended  tri-diagonal  matrix 
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0 
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(A5) 


where  the  M  x  M  block  matrices  ApB^andCj  are  functions  of  the  finite  difference  parameters  and  the 

Pedersen  conductivity  samples  Zpij .  The  matrices  Ai  and  Di  are  needed  for  periodic  boundaries  in  the  y- 

direction.  The  string  of  matrices  [Di,  D2,  D3,  Dn_i,  DN]  allow  inclusion  of  an  additional  condition  on  the 
potential  such  as 


JJ  ®(x,  y)  dx  dy  =  0  . 


(A6) 


This  condition  arises  when  the  addition  of  a  constant  to  a  solution  also  yields  solution.  The  nonuniqueness 
occurs  if  the  boundary  conditions  are  completely  periodic  and/or  specified  by  constant  derivatives  (i.e., 
Neumann).  With  these  types  of  boundary  conditions,  (A6)  prevents  the  square  matrix  in  (A5)  from  being 
singular  and  a  numerical  solution  can  be  obtained.  Finally,  the  right  side  of  (Al)  and  boundary  conditions  are 
contained  in  the  linear  arrays  Yj . 

The  algorithm  for  solving  (A5)  is  a  generalization  of  the  Thomas  Algorithm  for  scalar  tridiagonal  systems  [5]. 
Initialize  with  new  matrix  variables 


oq  =  Bj1,  Sx  =  oq  •  Y19  a2  =  cq  C19  bx  =  -a,  •  Ax  (A7) 

Continue  with  the  equations 

=  aM  •  CM,  oq  =  [Bj  -  Ai  •  aj"1,  =  cq  •  [Yj  -  A-S^J,  =  -oq  •  AM  •  b{_t  (A8) 

for  the  index  in  the  range  i  =  2, ... ,  N  - 1 .  The  next  operations  define  a  new  set  of  variables  starting  with 
S^=0,  b^  =  I ,  where  I  is  the  M  x  M  identity  matrix,  and  continue  with 


Sn_.  =  sN_, 


aN-i+l^N-i+l’  ^N  i  ^N  i 


aN-i+l^N-i+l 


(A9) 
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for  the  index  “i”  in  the  range  i  =  l,...,N  —  1.  The  solution  for  j  =  N  is  given  by 


(A10) 


i=l  J  L  i=l 


The  remaining  solution  vectors  are  found  from 


x.  =s;  +b;  -xN 


(All) 


where  the  index  “i”  is  given  by  the  range  i  =  1, ...  ,N  —  1.  This  algorithm  was  used  to  provide  the  numerical 
potential  solutions  illustrated  in  Figure  5  and  the  data  given  in  Tables  II  and  III. 
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